---
title: "Empirical evidence of mental health risks posed by climate change"
subtitle: Supplementary information replication code
author: Nick Obradovich, Robyn Migliorini, Martin P. Paulus, and Iyad Rahwan
output: pdf_document
---

<!--
Rscript -e "rmarkdown::render('./si_replication_code.Rmd')"
-->

```{r opts, echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
knitr::opts_chunk$set(fig.path='.', warning=F, message=F, fig.retina=NULL, cache=T, cache.lazy=F, autodep=T, echo=F, eval=T)
```

```{r packages}
# packages----
library(bit64) # long integers
library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
library(knitr) #knitting
library(pander) #pandering
library(rasterVis) # plotting rasters
library(GGally) # for correlation plots
library(geosphere) # spherical distance
library(Rcpp) # for C++
library(RcppArmadillo) # for C++ (and Conley errors)
library(matrixStats) # matrix calculations
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(viridis) # for colors
library(Cairo) # pdf device
library(gtable) # for gtable filter
library(devtools)
#devtools::install_github("wilkelab/cowplot")
library(cowplot)
library(captioner)

# set cores for parallel processing----
registerDoMC(20)
```

```{r functions}
# functions for grabbing p-values and coefficients from summaries of felm models----
coefs <- function(reg, name) {
    if (round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3) == 0) {
        format(round(reg$coefficients[rownames(reg$coefficients)==name, 1], 4), scientific=FALSE) 
    }
    else round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3)
}
pval <- function(reg, name) {
  if (round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3) >=0.001) 
  {round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3)} 
  else paste("< 0.001")
}

# captioner setting----
fig_num <- captioner(prefix="Fig.")
tab_num <- captioner(prefix="Table")
eq_num <- captioner(prefix="Equation")
sfig_num <- captioner(prefix="Fig. S", auto_space=FALSE)
stab_num <- captioner(prefix="Table S", auto_space=FALSE)
seq_num <- captioner(prefix="Equation S", auto_space=FALSE)

# binned plot----
source('./binned_plot_function.R')
```

```{r data}
load('./climate_mental_health.Rdata')
```

<!-- figures-->

```{r fig.respondents, eval=FALSE}
# plot----
p <- ggplot(brfss[year(time) < 2012,]) + 
        geom_point(aes(x=time, y=name), alpha=0.01) +
        coord_cartesian(expand = FALSE) +
          theme(
                plot.title = element_text(hjust=0.025),
                axis.title.x = element_blank(),
                axis.title.y = element_blank(),
                axis.text.y = element_text(size=6, face="bold"),
                axis.ticks.y = element_blank(),
                axis.text.x = element_text(size=20, face="bold"),
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank(),
                panel.background = element_blank(),
                axis.line = element_line(),
                legend.title=element_blank()
          )

png(filename = './respondents.png', width=15, height=19, units = "in", type="cairo", res=150)
p
dev.off()
```

```{r respondent, out.width="100%", fig.cap = paste0(sfig_num("respondents", display="cite"),": Plot of all respondent-days across the 263 unique cities in the sample. This figure plots a point by city for each of the nearly two million respondents in our sample.")}
knitr::include_graphics("./respondents.png")
```

```{r}
respondents <- sfig_num(name = "respondents", caption = "")
```

```{r fig.decon, eval=FALSE}
# example unit to demonstrate effect of deconvolving fixed effects on temperature data----
de <- as.data.table(demeanlist(mtx = list(brfss$avg.tmax), fl = list(brfss$mmsa, brfss$stateyr, as.factor(brfss$time))))
setnames(de, c('tmax.demean.all'))
demeaned <- cbind(brfss, de)
de <- as.data.table(demeanlist(mtx = list(brfss$avg.tmax), fl = list(brfss$mmsa)))
setnames(de, c('tmax.demean.unitonly'))
demeaned <- cbind(demeaned, de)

# Albuquerque NM----
p1 <- ggplot(demeaned[name=="Albuquerque NM",], aes(x=time)) + 
    ylab("30-Day Avgerage Max. Temp. in \u2103") +
    xlab("Albuquerque, NM") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=avg.tmax, color='black'), size=0.5, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.5, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color="#fc8d62"), size=0.5, alpha=1) +
    scale_color_manual(limits=c("black","blue4","#fc8d62"), values=c("black","blue4","#fc8d62"), labels=c("No Fixed Effects", "City Fixed Effects", "City, Date, and State*Year Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.position = c(0.05,0.15),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )

# Denver CO----
p2 <- ggplot(demeaned[name=="Denver CO",], aes(x=time)) + 
    ylab("30-Day Avgerage Max. Temp. in \u2103") +
    xlab("Denver, CO") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=avg.tmax, color='black'), size=0.5, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.5, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color="#fc8d62"), size=0.5, alpha=1) +
    scale_color_manual(limits=c("black","blue4","#fc8d62"), values=c("black","blue4","#fc8d62"), labels=c("No Fixed Effects", "City Fixed Effects", "City, Date, and State*Year Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.position = c(0.05,0.15),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )

# Louisville KY----
p3 <- ggplot(demeaned[name=="Louisville KY",], aes(x=time)) + 
    ylab("30-Day Avgerage Max. Temp. in \u2103") +
    xlab("Louisville, KY") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=avg.tmax, color='black'), size=0.5, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.5, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color="#fc8d62"), size=0.5, alpha=1) +
    scale_color_manual(limits=c("black","blue4","#fc8d62"), values=c("black","blue4","#fc8d62"), labels=c("No Fixed Effects", "City Fixed Effects", "City, Date, and State*Year Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.position = c(0.05,0.15),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )

# Phoenix AZ----
p4 <- ggplot(demeaned[name=="Phoenix AZ",], aes(x=time)) + 
    ylab("30-Day Avgerage Max. Temp. in \u2103") +
    xlab("Phoenix, AZ") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=avg.tmax, color='black'), size=0.5, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.5, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color="#fc8d62"), size=0.5, alpha=1) +
    scale_color_manual(limits=c("black","blue4","#fc8d62"), values=c("black","blue4","#fc8d62"), labels=c("No Fixed Effects", "City Fixed Effects", "City, Date, and State*Year Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.position = c(0.05,0.15),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )

# justy----
test <- ggplot(demeaned[name=="Albuquerque NM",], aes(x=time)) + 
    ylab("30-Day Max. Temp. in \u2103") +
    xlab("Albuquerque, NM") +
    geom_hline(yintercept = 0, linetype="dashed") +
    geom_line(aes(y=avg.tmax, color='black'), size=0.5, alpha=0.5) + 
    geom_line(aes(y=tmax.demean.unitonly, color='blue4'), size=0.5, alpha=0.5) +
    geom_line(aes(y=tmax.demean.all, color="#fc8d62"), size=0.5, alpha=1) +
    scale_color_manual(limits=c("black","blue4","#fc8d62"), values=c("black","blue4","#fc8d62"), labels=c("No Fixed Effects", "City Fixed Effects", "City, Date, and State*Year Fixed Effects")) +
    guides(color = guide_legend(override.aes = list(size=5))) +
    coord_cartesian(ylim=c(-35, 45), expand=TRUE) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_text(size=25, angle=90, face="bold"),
          axis.text.y = element_text(size=20, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position ="none",
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
justy <- gtable_filter(ggplotGrob(test), 'axis-l|ylab', trim=F)
justy <- ggdraw(justy)

# plot----
png(filename = "./deconvolved_avg_tmax.png", width=18, height=10, units = "in", type="cairo", res=300)
grid.arrange(
    grobs = list(
        arrangeGrob(grobs=list(justy, arrangeGrob(p1, p2, ncol=2)), widths=c(0.05, 0.95), ncol=2, nrow=1),
        arrangeGrob(grobs=list(justy, arrangeGrob(p3, p4, ncol=2)), widths=c(0.05, 0.95), ncol=2, nrow=1)
    )
)
grid.text("a", x=unit(0.065, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.55, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("c", x=unit(0.065, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("d", x=unit(0.55, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r decon, out.width="100%", fig.cap=paste0(sfig_num("decon", display="cite"), ": Deconvolved thirty-day running average temperature variation. Each panel of this figure depicts three series. The first is the raw thirty-day running average temperature observations (No Fixed Effects), the second is the remaining thirty-day average temperature variation after deconvolving the raw data around only city-level fixed effects, the third is the remaining variation after deconvolving the raw data around our full set of city, date, and state-by-year fixed effects. As can be seen, trends are removed and the remaining identifying variation is stationary (not systematically trending over time) and centered around zero for each individual unit. Panel (a) plots this relationship for Albuquerque, NM. Panel (b) plots this relationship for Denver, CO. Panel (c) plots this relationship for Louisville, KY. Panel (d) plots this relationship for Phoenix, AZ.")}
knitr::include_graphics("./deconvolved_avg_tmax.png")
```

```{r}
decon <- sfig_num(name = "decon", caption = "")
```

```{r fig.bins, eval=FALSE}
# varying bins----
## 5 degree bins
brfss[, cuttmax := cut(avg.tmax, breaks=c(-Inf,seq(0, 30, by=5),Inf))] # create factor variable for regresison
brfss[, cuttmax := relevel(cuttmax, ref="(10,15]")] # set reference category
brfss[, cutprcp := cut(prcp.days, breaks=c(-Inf,seq(0, 25, by=5),Inf))] # create factor variable for regresison
brfss[, cutprcp := relevel(cutprcp, ref="(-Inf,0]")] # set reference category
cut5 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)

## 2 degree bins
brfss[, cuttmax := cut(avg.tmax, breaks=c(-Inf,seq(0, 30, by=2),Inf))] # create factor variable for regresison
brfss[, cuttmax := relevel(cuttmax, ref="(12,14]")] # set reference category
brfss[, cutprcp := cut(prcp.days, breaks=c(-Inf,seq(0, 24, by=2),Inf))] # create factor variable for regresison
brfss[, cutprcp := relevel(cutprcp, ref="(-Inf,0]")] # set reference category
cut2 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)

## 1 degree bins
brfss[, cuttmax := cut(avg.tmax, breaks=c(-Inf,seq(0, 30, by=1),Inf))] # create factor variable for regresison
brfss[, cuttmax := relevel(cuttmax, ref="(12,13]")] # set reference category
brfss[, cutprcp := cut(prcp.days, breaks=c(-Inf,seq(0, 25, by=1),Inf))] # create factor variable for regresison
brfss[, cutprcp := relevel(cutprcp, ref="(-Inf,0]")] # set reference category
cut1 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE,
          psdef = FALSE)

## set back to original bins
brfss[, cuttmax := cut(avg.tmax, breaks=c(-Inf,seq(0, 30, by=5),Inf))] # create factor variable for regresison
brfss[, cuttmax := relevel(cuttmax, ref="(10,15]")] # set reference category
brfss[, cutprcp := cut(prcp.days, breaks=c(-Inf,seq(0, 25, by=5),Inf))] # create factor variable for regresison
brfss[, cutprcp := relevel(cutprcp, ref="(-Inf,0]")] # set reference category

# tmax plots----
cut5.plot <- bin.plot(felm.est = cut5,
                     plotvar = "tmax",
                     breaks = 5,
                     omit = c(10,15),
                     ylimit = c(-1, 3),
                     xlabel = "30-Day Avg. Max. Temp. in \u2103",
                     ylabel = "P(Mental Health Issue)\n(Percentage Points)",
                     y.axis.title = element_blank(),
                     y.axis.text = element_blank(),
                     y.axis.line = element_blank(),
                     y.axis.ticks = element_blank(),
                     panel = ""
                     )
cut5.plot <- cut5.plot + theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, angle=0, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=8, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank()
                                )
hist <- axis_canvas(cut5.plot, axis = "x") + 
        geom_histogram(data = prism, aes(x = avg.tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
cut5.plot <- insert_xaxis_grob(cut5.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
cut5.plot <- ggdraw(cut5.plot)


cut2.plot <- bin.plot(felm.est = cut2,
                     plotvar = "tmax",
                     breaks = 2,
                     omit = c(12,14),
                     ylimit = c(-1, 3),
                     xlabel = "30-Day Avg. Max. Temp. in \u2103",
                     ylabel = "P(Mental Health Issue)\n(Percentage Points)",
                     y.axis.title = element_blank(),
                     y.axis.text = element_blank(),
                     y.axis.line = element_blank(),
                     y.axis.ticks = element_blank(),
                     panel = ""
                     )
cut2.plot <- cut2.plot + theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, angle=0, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=8, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank()
                                )
hist <- axis_canvas(cut2.plot, axis = "x") + 
        geom_histogram(data = prism, aes(x = avg.tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
cut2.plot <- insert_xaxis_grob(cut2.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
cut2.plot <- ggdraw(cut2.plot)


cut1.plot <- bin.plot(felm.est = cut1,
                     plotvar = "tmax",
                     breaks = 1,
                     omit = c(12,13),
                     ylimit = c(-1, 3),
                     xlabel = "30-Day Avg. Max. Temp. in \u2103",
                     ylabel = "P(Mental Health Issue)\n(Percentage Points)",
                     y.axis.title = element_blank(),
                     y.axis.text = element_blank(),
                     y.axis.line = element_blank(),
                     y.axis.ticks = element_blank(),
                     panel = ""
                     )
cut1.plot <- cut1.plot + theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, angle=0, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=8, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank()
                                )
hist <- axis_canvas(cut1.plot, axis = "x") + 
        geom_histogram(data = prism, aes(x = avg.tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
cut1.plot <- insert_xaxis_grob(cut1.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
cut1.plot <- ggdraw(cut1.plot)

## y-axis
ylabel <- bin.plot(felm.est = cut5,
                     plotvar = "tmax",
                     breaks = 5,
                     omit = c(10,15),
                     ylimit = c(-1, 3),
                     xlabel = "30-Day Avg. Max. Temp. in \u2103",
                     ylabel = "P(Mental Health Issue)\n(Percentage Points)",
                     # y.axis.title = element_blank(),
                     # y.axis.text = element_blank(),
                     # y.axis.line = element_blank(),
                     # y.axis.ticks = element_blank(),
                     panel = ""
                     )
hist <- axis_canvas(ylabel, axis = "x") + 
        geom_histogram(data = prism, aes(x = tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
ylabel <- insert_xaxis_grob(ylabel, hist, position = "bottom", height = grid::unit(0.05, "null"))
ylabel <- gtable_filter(ylabel, 'axis-l|ylab', trim=F)
ylabel <- ggdraw(ylabel)

# prcp plots----
## 5 degree bins
cutp5.plot <- bin.plot(felm.est = cut5,
                     plotvar = "prcp",
                     breaks = 5,
                     omit = c(0,0),
                     ylimit = c(-1, 3),
                     xlabel = "Monthly Days of Precipitation",
                     ylabel = "P(Mental Health Issue)\n(Percentage Points)",
                     y.axis.title = element_blank(),
                     y.axis.text = element_blank(),
                     y.axis.line = element_blank(),
                     y.axis.ticks = element_blank(),
                     panel = ""
                     )
cutp5.plot <- cutp5.plot + theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, angle=0, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=8, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank()
                                )
hist <- axis_canvas(cutp5.plot, axis = "x") + 
        geom_histogram(data = prism, aes(x = avg.tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
cutp5.plot <- insert_xaxis_grob(cutp5.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
cutp5.plot <- ggdraw(cutp5.plot)

## 2 degree bins
cutp2.plot <- bin.plot(felm.est = cut2,
                     plotvar = "prcp",
                     breaks = 2,
                     omit = c(0,0),
                     ylimit = c(-1, 3),
                     xlabel = "Monthly Days of Precipitation",
                     ylabel = "P(Mental Health Issue)\n(Percentage Points)",
                     y.axis.title = element_blank(),
                     y.axis.text = element_blank(),
                     y.axis.line = element_blank(),
                     y.axis.ticks = element_blank(),
                     panel = ""
                     )
cutp2.plot <- cutp2.plot + theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, angle=0, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=8, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank()
                                )
hist <- axis_canvas(cutp2.plot, axis = "x") + 
        geom_histogram(data = prism, aes(x = avg.tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
cutp2.plot <- insert_xaxis_grob(cutp2.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
cutp2.plot <- ggdraw(cutp2.plot)

## 1 degree bins
cutp1.plot <- bin.plot(felm.est = cut1,
                     plotvar = "prcp",
                     breaks = 1,
                     omit = c(0,0),
                     ylimit = c(-1, 3),
                     xlabel = "Monthly Days of Precipitation",
                     ylabel = "P(Mental Health Issue)\n(Percentage Points)",
                     y.axis.title = element_blank(),
                     y.axis.text = element_blank(),
                     y.axis.line = element_blank(),
                     y.axis.ticks = element_blank(),
                     panel = ""
                     )
cutp1.plot <- cutp1.plot + theme(plot.title = element_text(face="bold", size=40),
                              axis.title.x = element_text(size=25, angle=0, face="bold"),
                              axis.title.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.text.x = element_text(size=8, face="bold"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(),
                              axis.line.x = element_line(),
                              axis.line.y = element_blank(),
                              axis.ticks.y = element_blank(),
                              legend.title=element_blank()
                                )
hist <- axis_canvas(cutp1.plot, axis = "x") + 
        geom_histogram(data = prism, aes(x = avg.tmax), bins=30, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
cutp1.plot <- insert_xaxis_grob(cutp1.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
cutp1.plot <- ggdraw(cutp1.plot)

# plot----
png(filename = './varying_bins.png', width=24, height=12, units = "in", type="cairo", res=300)
grid.arrange(
    arrangeGrob(ylabel, arrangeGrob(cut1.plot, cut2.plot, cut5.plot, ncol=3), ncol=2, widths=c(0.06, 0.94)),
    arrangeGrob(ylabel, arrangeGrob(cutp1.plot, cutp2.plot, cutp5.plot, ncol=3), ncol=2, widths=c(0.06, 0.94)),
    nrow=2
)
grid.text("a", x=unit(0.07, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.39, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("c", x=unit(0.70, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("d", x=unit(0.07, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("e", x=unit(0.39, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("f", x=unit(0.70, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=40,font=2))
dev.off()
```

```{r bins, out.width="100%", fig.cap=paste0(sfig_num("binsizes", display="cite"), ": Varying bin sizes of average thirty-day maximum temperatures and number of precipitation days. This figure replicates the models from Figure 1 in the main text with varying bin sizes across 1, 2, and 5 unit bins. Panel (a) plots the estimates from the main text model with 1 degree Celsius temperature bins. Panel (b) uses 2 degree Celsius temperature bins. Panel (c) uses 5 degree Celsius bins. Panel (d) plots the estimates from the main text model with 1-day precipitation day bins. Panel (e) uses 2-day bins. Panel (f) uses 5-day bins. Results are robust across choice of temperature and precipitation bin size.")}
knitr::include_graphics("./varying_bins.png")
```

```{r}
binsizes <- sfig_num(name = "binsizes", caption = "")
```

```{r fig.par, eval=FALSE}
# parametric regression----
brfss[, avg.tmax2 := avg.tmax^2]
brfss[, prcp.days2 := prcp.days^2]
quad <- felm(any_ment ~ avg.tmax + avg.tmax2 + prcp.days + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)

# temp----
t <- brfss[, seq(from=-10, to=40, by=0.001)]
d.t <- quad$coef[1]*t + quad$coef[2]*t^2
t.data <- as.data.table(cbind(t, d.t))

t.plot <- ggplot(t.data) + 
    ylab("P(Report Mental Health Issue)\n(Percentage Points)") +
    xlab("30-Day Avg. Max. Temp. in \u2103") +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=d.t), color="#fc8d62", size=1.5) + 
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    coord_cartesian(xlim=c(-10,40)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, face="bold"),
          axis.title.y = element_text(size=25, face="bold"),
          axis.text.y = element_text(size=20, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_line(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position = "none"
    )
hist <- axis_canvas(t.plot, axis = "x") + 
        geom_histogram(data = brfss, aes(x = avg.tmax), bins=20, fill='#8da0cb', color='gray60', alpha=0.5, size=0.1)
t.plot <- insert_xaxis_grob(t.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
t.plot <- ggdraw(t.plot)

# margn effect temp----
# get cluster robust variance covariance matrices
quad.vcv <- quad$clustervcv

# get coefficients
quad.tcoef <- quad$coef[1]
quad.t2coef <- quad$coef[2]

# calculate marginal effects
marg <- data.table(tmax = seq(from=-10, to=40, by=.001))
marg[, quad.delta := quad.tcoef + 2*quad.t2coef*tmax]

# compute variance of marginal effect
marg[, quad.se := sqrt(quad.vcv["avg.tmax","avg.tmax"] + 4*(tmax^2)*quad.vcv["avg.tmax2","avg.tmax2"] + 4*tmax*quad.vcv[1,2])]

# 95% confidence bounds
marg[, quad.upr := quad.delta + 1.96*quad.se]
marg[, quad.lwr := quad.delta - 1.96*quad.se]

quad.plot <- ggplot(data=marg) + 
                ylab("Estimated Marginal Effect of +1\u2103 on\n P(Report Mental Health Issue) in Percentage Points") +
                xlab("30-Day Avg. Max. Temp. in \u2103") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_ribbon(aes(x=tmax, ymin=quad.lwr, ymax=quad.upr), fill="#fc8d62", alpha=0.25) + 
                geom_line(aes(x=tmax, y=quad.delta), color="#fc8d62", size=1.5) +
                scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )
hist <- axis_canvas(quad.plot, axis = "x") + 
        geom_histogram(data = brfss, aes(x = avg.tmax), bins=20, fill='#8da0cb', color='gray60', alpha=0.5, size=0.1)
quad.plot <- insert_xaxis_grob(quad.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
quad.plot <- ggdraw(quad.plot)

# prcp----
p <- brfss[, seq(from=min(prcp.days, na.rm=T), to=max(prcp.days, na.rm=T), by=0.001)]
d.p <- quad$coef[3]*p
p.data <- as.data.table(cbind(p, d.p))
p.data[, upr := d.p + quad$cse[3]*1.96]
p.data[, lwr := d.p - quad$cse[3]*1.96]

p.plot <- ggplot(p.data) + 
    ylab("P(Report Mental Health Issue)\n(Percentage Points)") +
    xlab("Monthly Days of Precipitation") +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_ribbon(aes(x=p, ymin=lwr, ymax=upr), fill="#8da0cb", alpha=0.25) + 
    geom_line(aes(x=p, y=d.p), color="#8da0cb", size=1.5) + 
    scale_x_continuous(breaks=c(0, 10, 20, 30)) +
    coord_cartesian(xlim=c(0,30)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, face="bold"),
          axis.title.y = element_text(size=25, face="bold"),
          axis.text.y = element_text(size=20, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_line(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position = "none"
    )
hist <- axis_canvas(p.plot, axis = "x") + 
        geom_histogram(data = brfss, aes(x = prcp.days), bins=20, fill='#8da0cb', color='gray60', alpha=0.5, size=0.1)
p.plot <- insert_xaxis_grob(p.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
p.plot <- ggdraw(p.plot)

# plot----
png(filename = './parametric_regs.png', width=22, height=10, units = "in", type="cairo", res=300)
grid.arrange(t.plot, quad.plot, p.plot, ncol=3)
grid.text("a", x=unit(0.075, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.415, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("c", x=unit(0.74, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
dev.off()
```

```{r parametrics, out.width="100%", fig.cap=paste0(sfig_num("parametric", display="cite"), ": Parametrically modeling maximum temperature and precipitation also uncovers a positive effect on mental health issues. This figure draws on the estimation of the linear probability model in Equation 1 and plots the results of parametrically modeling the relationships between our meteorological variables and our mental health outcome variable. Panel (a) displays the quadratic relationship between average maximum temperatures over the respective thirty-day window and reports of mental health issues. Panel (b) plots the marginal effects and 95\\% confidence intervals associated with the fitted values from panel (a). Panel (c) depicts the linear relationship between number of days with precipitation and probability of reporting mental health issues. The probability of reporting mental health issues increases linearly with precipitation days and quadratically with higher average maximum temperatures over the thirty-day window prior to participant response. Shaded error regions represent 95\\% confidence intervals.")}
knitr::include_graphics("./parametric_regs.png")
```

```{r}
parametric <- sfig_num(name = "parametric", caption = "")
```

```{r fig.daybins, eval=FALSE}
# reg----
day <- felm(any_ment ~ cutprcp + tmaxneginf + tmax5 + tmax10 + tmax20 + tmax25 + tmax30 + tmax35 + tmaxposinf 
            + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)

# plot----
day.plot <- bin.plot(felm.est = day,
                     plotvar = "tmax",
                     breaks = 5,
                     omit = c(10,15),
                     ylimit = c(-0.05, 0.11),
                     xlabel = "Count of Daily Max. Temp. in \u2103",
                     ylabel = "P(Report Mental Health Issue)\n(Percentage Points)",
                     panel = ""
                     )
hist <- axis_canvas(day.plot, axis = "x") + 
        geom_histogram(data = prism, aes(x = tmax), bins=20, fill="#8da0cb", color='gray60', alpha=0.5, size=0.1)
day.plot <- insert_xaxis_grob(day.plot, hist, position = "bottom", height = grid::unit(0.05, "null"))
day.plot <- ggdraw(day.plot)

# plots----
png(filename = './daybin_plot.png', width=8, height=8, units = "in", type="cairo", res=300)
day.plot
dev.off()
```

```{r daybins, out.width="100%", fig.cap=paste0(sfig_num("daybinfigs", display="cite"), ": Acconting for the distribution of daily maximum temperatures also shows positive effect of hotter temperatures on reports of mental health issues. This figure draws on the estimation of the linear probability model in Equation 1. It pulls in the full distribution of daily maximum temperatures over the thirty-day window and plots the marginal effects of an added day in each daily temperature bin on the probability of reporting a mental health issue. Coefficient estimates represent the effect on probability of reporting mental health issues of having one additional day out of the thirty-day window falling in the specified maximum temperature bin. The probability of mental health issues steadily increases with each additional day that exceeds 10-15$^\\circ{}$C (50-60F).")}
knitr::include_graphics("./daybin_plot.png")
```

```{r}
daybinfigs <- sfig_num(name = "daybinfigs", caption = "")
```

```{r fig.loo, eval=FALSE}
# leave one out regression----
# a <- foreach (i = brfss[, unique(state)]) %do% {
#     a <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
#                  | time + mmsa + stateyr
#                  | 0
#                  | state
#                  , data=brfss[state!=i],
#                  exactDOF=TRUE))
#     a <- data.table(coef25=a$coefficients[12, 1], coef30=a$coefficients[13, 1], cse25=a$coefficients[12, 2], cse30=a$coefficients[13, 2], state=i)
#     return(a)
# }
# a <- rbindlist(a)
# saveRDS(a, "./leave_one_out_coefficients.rds")
a <- readRDS("./leave_one_out_coefficients.rds")
a[, upr30 := coef30 + 1.96*cse30]
a[, lwr30 := coef30 - 1.96*cse30]
a[, upr25 := coef25 + 1.96*cse25]
a[, lwr25 := coef25 - 1.96*cse25]

# coefficient plot----
loo.plot25 <- ggplot(data=a) + 
                ylab("Effect of 25-30\u2103 on\n P(Report Mental Health Issue)") +
                xlab("Omitted State") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_linerange(aes(ymax=upr25, ymin=lwr25, x=state), color="#8da0cb", size=1.25) +
                geom_point(aes(x=state, y=coef25), color="#fc8d62", size=3) +
                coord_cartesian(ylim=c(0,2)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_blank(),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_blank(),
                      axis.ticks.y = element_blank(),
                      axis.ticks.x = element_blank(),
                      axis.line.x = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

loo.plot30 <- ggplot(data=a) + 
                ylab("Effect of >30\u2103 on\n P(Report Mental Health Issue)") +
                xlab("Omitted State") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_linerange(aes(ymax=upr30, ymin=lwr30, x=state), color="#8da0cb", size=1.25) +
                geom_point(aes(x=state, y=coef30), color="#fc8d62", size=3) +
                coord_cartesian(ylim=c(0,2)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

# plot----
png(filename = './leave_one_out.png', width=23, height=12, units = "in", type="cairo", res=300)
grid.arrange(loo.plot25, loo.plot30, nrow=2)
grid.text("a", x=unit(0.07, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.07, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=40,font=2))
dev.off()
```

```{r loos, out.width="100%", fig.cap=paste0(sfig_num("loo", display="cite"), ": Results are robust to omission of observations from particular states. This figure plots the results of the estimation of the linear probability model in Equation 1 across fifty iterated regressions, each dropping one state per each estimation. The goal of this estimation is to examine whether any one particular state is substantially driving the results we report in the main text. As can be seen across both panel (a) (displays the coefficient on the 25-30C bin) and panel (b) (displays the coefficient on the >30C bin), our parameter estimates are quite consistent, regardless of which state's respondents are omitted from the analysis. Shaded error regions represent 95\\% confidence intervals.")}
knitr::include_graphics("./leave_one_out.png")
```

```{r}
loo <- sfig_num(name = "loo", caption = "")
```

```{r fig.div.loo, eval=FALSE}
# leave one out regression----
# a <- foreach (i = brfss[, unique(division)]) %do% {
#     a <- summary(felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
#                  | time + mmsa + stateyr
#                  | 0
#                  | state
#                  , data=brfss[division!=i],
#                  exactDOF=TRUE))
#     a <- data.table(coef25=a$coefficients[12, 1], coef30=a$coefficients[13, 1], cse25=a$coefficients[12, 2], cse30=a$coefficients[13, 2], division=i)
#     return(a)
# }
# a <- rbindlist(a)
# saveRDS(a, "./division_leave_one_out_coefficients.rds")
a <- readRDS("./division_leave_one_out_coefficients.rds")
a[, upr30 := coef30 + 1.96*cse30]
a[, lwr30 := coef30 - 1.96*cse30]
a[, upr25 := coef25 + 1.96*cse25]
a[, lwr25 := coef25 - 1.96*cse25]

# coefficient plot----
loo.plot25 <- ggplot(data=a) + 
                ylab("Effect of 25-30\u2103 on\n P(Report Mental Health Issue)") +
                xlab("Omitted Census Division") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_linerange(aes(ymax=upr25, ymin=lwr25, x=division), color="#8da0cb", size=1.25) +
                geom_point(aes(x=division, y=coef25), color="#fc8d62", size=3) +
                coord_cartesian(ylim=c(0,2)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_blank(),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_blank(),
                      axis.ticks.y = element_blank(),
                      axis.ticks.x = element_blank(),
                      axis.line.x = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

loo.plot30 <- ggplot(data=a) + 
                ylab("Effect of >30\u2103 on\n P(Report Mental Health Issue)") +
                xlab("Omitted Census Division") +
                geom_hline(aes(yintercept=0), linetype="dashed") +
                geom_linerange(aes(ymax=upr30, ymin=lwr30, x=division), color="#8da0cb", size=1.25) +
                geom_point(aes(x=division, y=coef30), color="#fc8d62", size=3) +
                coord_cartesian(ylim=c(0,2)) +
                theme(plot.title = element_text(face="bold", size=40),
                      axis.title.x = element_text(size=25, angle=0, face="bold"),
                      axis.title.y = element_text(size=25, face="bold"),
                      axis.text.y = element_text(size=20, face="bold"),
                      axis.text.x = element_text(size=20, face="bold"),
                      axis.ticks.y = element_blank(),
                      panel.grid.major = element_blank(),
                      panel.grid.minor = element_blank(),
                      panel.background = element_blank(),
                      legend.title=element_blank(),
                      legend.position = "none"
                )

# plot----
png(filename = './division_leave_one_out.png', width=25, height=12, units = "in", type="cairo", res=300)
grid.arrange(loo.plot25, loo.plot30, nrow=2)
grid.text("a", x=unit(0.07, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.07, "npc"), y=unit(0.47, "npc"), gp=gpar(fontsize=40,font=2))
dev.off()
```

```{r div.loos, out.width="100%", fig.cap=paste0(sfig_num("div.loo", display="cite"), ": Results are robust to omission of observations from particular US Census divisions. This figure plots the results of the estimation of the linear probability model in Equation 1 across nine iterated regressions, each dropping one census division per each estimation. The goal of this estimation is to examine whether any one particular division is substantially driving the results we report in the main text. As can be seen across both panel (a) (displays the coefficient on the 25-30C bin) and panel (b) (displays the coefficient on the >30C bin), our parameter estimates are moderately consistent, regardless of which division's respondents are omitted from the analysis. Omission of the Pacific census division attenuates effect sizes, indicating cities in the Pacific division have more pronounced responses to these hotter temperatures than average. Shaded error regions represent 95\\% confidence intervals.")}
knitr::include_graphics("./division_leave_one_out.png")
```

```{r}
div.loo <- sfig_num(name = "div.loo", caption = "")
```

```{r fig.ldperiods, eval=FALSE}
# plot of all the different combinations of long-difference periods----
ldcomb[, post.min := as.numeric(substr(post, start=1, stop=4))]
ldcomb[, post.max := as.numeric(substr(post, start=nchar(post)-3, stop=nchar(post)))]
ldcomb[, pre.min := as.numeric(substr(pre, start=1, stop=4))]
ldcomb[, pre.max := as.numeric(substr(pre, start=nchar(pre)-3, stop=nchar(pre)))]
ldcomb[, iter := .I]

p <- ggplot() + 
    geom_linerange(data=ldcomb, aes(x=iter, ymin=pre.min, ymax=pre.max + 1), color="#fc8d62", size=0.5) +
    geom_linerange(data=ldcomb, aes(x=iter, ymin=post.min, ymax=post.max + 1), color="#8da0cb", size=0.5) +
    scale_x_continuous(breaks=seq(from=0, to=250, by=50)) +
    scale_y_continuous(breaks=c(2002:2012)) +
    xlab("Iteration") +
    ylab("Date Range") +
    coord_flip() +
    theme(
        axis.title.x = element_text(size=25, face="bold"),
        axis.title.y = element_text(size=25, face="bold"),
        axis.text.y = element_text(size=20, face="bold"),
        axis.text.x = element_text(size=20, face="bold"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(),
        legend.title=element_blank()
    )

png(filename = './ldperiods.png', width=15, height=8, units = "in", type="cairo", res=300)
p
dev.off()
```

```{r ldperiods, out.width="100%", fig.cap=paste0(sfig_num("ldperiods", display="cite"), ": Plot of all \\textit{a} periods and \\textit{b} periods in the long-differences permutation exercise. This figure plots a line-range for each \\textit{a} period and \\textit{b} period for each iteration of the permutation of the long-differences model. \\textit{a} periods are plotted in orange while \\textit{b} periods are plotted in blue. \\textit{a} period years are strictly five or more years prior to \\textit{b} period years.")}
knitr::include_graphics("./ldperiods.png")
```

```{r}
ldperiods <- sfig_num(name = "ldperiods", caption = "")
```

```{r fig.menthlth, eval=FALSE}
# plot----
p <- ggplot(brfss) + 
        geom_histogram(aes(x=menthlth, y=..density..), bins=30, color="blue4", fill="#fc8d62") +
        scale_x_continuous(breaks=seq(0, 30, 5)) +
        xlab("Days with Mental Health Difficulties") +
        ylab("Density") +
        theme(
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_text(size=25, angle=90, face="bold"),
          axis.text.y = element_text(size=20, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position ="none",
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
        )
png(filename = './menthlth.png', width=10, height=5, units = "in", type="cairo", res=300)
p
dev.off()
```

```{r menthlth, out.width="100%", fig.cap = paste0(sfig_num("menthlths", display="cite"),": Histogram of number of days of mental health difficulties. This figure plots the histogram of the count variable measured by the BRFSS. As can be seen, the number of reported days are likely measured with error, as there is notable clustering around week and 5-day intervals.")}
knitr::include_graphics("./menthlth.png")
```

```{r}
menthlths <- sfig_num(name = "menthlths", caption = "")
```

```{r fig.warm, eval=FALSE}
# create distribution of projection data----
ft <- data.table(f2010 = values(bcsd.tmax[[1]]), f2050 = values(bcsd.tmax[[2]]), f2099 = values(bcsd.tmax[[3]]))
ft <- melt(ft, variable = "year", value="tcount")
ft[, year := as.integer(gsub("f", "", as.character(year)))]
fp <- data.table(f2010 = values(bcsd.prcp.days[[1]]), f2050 = values(bcsd.prcp.days[[2]]), f2099 = values(bcsd.prcp.days[[3]]))
fp <- melt(fp, variable = "year", value="pcount")
fp[, year := as.integer(gsub("f", "", as.character(year)))]

p1 <- ggplot(data=ft[!is.na(tcount),], aes(x=tcount, color=factor(year, levels=c("2010","2050","2099")))) +
        geom_line(stat="density", size=2) +
        scale_color_manual(values=c("gray50","orange","red4"),labels=c("2010","2050","2099")) +
      ylab("Density") +
      xlab("Days with Max. Temp >30\u2103") +
        guides(fill = guide_legend(override.aes = list(colour = NULL))) +
        theme(axis.title.x = element_text(size=25, angle=0, vjust=-0, face="bold"),
            axis.title.y = element_text(size=25,angle=90, vjust=1, face="bold"),
            axis.text.y = element_text(size=20,angle=0, vjust=.5, face="bold") ,
            axis.text.x = element_text(size=20, face="bold"),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            panel.background = element_blank(),
            axis.line = element_line(),
            legend.title=element_blank(),
            legend.position = c(.30, 0.35),
            legend.text=element_text(size=18, face="bold"),
            legend.key = element_rect(colour = "black")
        )

p2 <- ggplot(data=fp[!is.na(pcount),], aes(x=pcount, color=factor(year, levels=c("2010","2050","2099")))) +
        geom_line(stat="density", size=2) +
        scale_color_manual(values=c("gray50","steelblue","blue4"),labels=c("2010","2050","2099")) +
      ylab("Density") +
      xlab("Days with Measurable Precipitation") +
        guides(fill = guide_legend(override.aes = list(colour = NULL))) +
        theme(axis.title.x = element_text(size=25, angle=0, vjust=-0, face="bold"),
            axis.title.y = element_text(size=25,angle=90, vjust=1, face="bold"),
            axis.text.y = element_text(size=20,angle=0, vjust=.5, face="bold") ,
            axis.text.x = element_text(size=20, face="bold"),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            panel.background = element_blank(),
            axis.line = element_line(),
            legend.title=element_blank(),
            legend.position = c(.5, 0.15),
            legend.text=element_text(size=18, face="bold"),
            legend.key = element_rect(colour = "black")
        )

# plot----
png(filename = "./climate_distributions.png", width=15, height=7.5, units = "in", type="cairo", res=300)
grid.arrange(p1, p2, ncol=2)
grid.text("a", x=unit(0.11, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.6, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r warm, out.width="100%", fig.cap=paste0(sfig_num("warms", display="cite"), ": Projected increase in number of hot days and change in days with measurable precipitation in the US. This figure presents the 25km x 25km grid cell projections of warming and future precipitation patterns in the United States over this century along the high emissions RCP8.5 scenario. We use downscaled climatic model data from NASA's NEX and calculate the number of days that exceed 30$^\\circ{}$C as well as those days that have measurable precipitation (>1mm) across the 21 models in the ensemble. We take the ensemble average for 2010, 2050, and 2099 for each of these metrics across each grid-cell in the data. As can be seen, the days with measurable precipitation projected relationship appears mostly unchanged over the coming century. However, days with hot temperatures are projected to increase substantially by 2099.")}
knitr::include_graphics("./climate_distributions.png")
```

```{r}
warms <- sfig_num(name = "warms", caption = "")
```

```{r fig.climchange, eval=FALSE}
# change in future hot days----
## trim forecast rasters to remove outer NA values
year2010 <- trim(bcsd.tmax[[1]], padding=0, values=NA)
year2050 <- trim(bcsd.tmax[[2]], padding=0, values=NA)
year2099 <- trim(bcsd.tmax[[3]], padding=0, values=NA)

## get differences
diff2050 <- year2050 - year2010
diff2099 <- year2099 - year2010

## stack these together
forecast <- stack(diff2050, diff2099)
id <- c('2050','2099') # set year id
forecast <- setZ(forecast, id) # set year as z variable
names(forecast) <- c('2050','2099') # set names

## plot raster forecasts
theme = rasterTheme(region=c('gray95','gray60','#fc8d62','red4')) # define color palette

p <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               col=NA,
               margin=TRUE, # remove margin padding
               scales=list(draw=FALSE),
               layout=c(2,1), # horizontal layout
               names.attr=c("", ""), # inset titles
               colorkey=list(space="right", labels=list(cex=1.25), height=0.75), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE,
               par.strip.text=list(cex=2, lines=1.5),
               sub=list("\u0394 Days with Maxmimum Temperature >30°C", cex=2.25)
               
    )
p <- update(p, par.settings = 
                list(axis.line = list(col = 0), 
                     panel.background=list(col="transparent")
                     )
            ) # remove plot borders and add subtitle below legend
p <- p + layer(sp.lines(states, lwd=0.5, col="gray60")) # add us shapefiles

# future prcp days----
## trim forecast rasters to remove outer NA values
year2010 <- trim(bcsd.prcp.days[[1]], padding=0, values=NA)
year2050 <- trim(bcsd.prcp.days[[2]], padding=0, values=NA)
year2099 <- trim(bcsd.prcp.days[[3]], padding=0, values=NA)

## get differences
diff2050 <- year2050 - year2010
diff2099 <- year2099 - year2010

## stack these together
forecast <- stack(diff2050, diff2099)
id <- c('2050','2099') # set year id
forecast <- setZ(forecast, id) # set year as z variable
names(forecast) <- c('2050','2099') # set names

## plot raster forecasts
theme = rasterTheme(region=c('gray95','gray60','#8da0cb','blue4')) # define color palette

p1 <- levelplot(forecast, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               col=NA,
               margin=TRUE, # remove margin padding
               scales=list(draw=FALSE),
               layout=c(2,1), # horizontal layout
               names.attr=c("", ""), # inset titles
               colorkey=list(space="right", labels=list(cex=1.25), height=0.75), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE,
               par.strip.text=list(cex=2, lines=1.5),
               sub=list("\u0394 Days with Measurable Precipitation", cex=2.25)
               
    )
p1 <- update(p1, par.settings = 
                list(axis.line = list(col = 0), 
                     panel.background=list(col="transparent")
                     )
            ) # remove plot borders and add subtitle below legend
p1 <- p1 + layer(sp.lines(states, lwd=0.5, col="gray60")) # add us shapefiles

# plot----
png(filename = './climate_changes.png', width=12, height=8, units = "in", type="cairo", res=300)
grid.arrange( 
             p,
             p1,
             nrow=2, ncol=1, heights=c(0.5, 0.5)
             )
grid.text("a", x=unit(0.04, "npc"), y=unit(0.94, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("b", x=unit(0.04, "npc"), y=unit(0.44, "npc"), gp=gpar(fontsize=40,font=2))
grid.text("\u0394 2050", x=unit(0.25, "npc"), y=unit(0.75, "npc"), gp=gpar(fontsize=25,font=2, col="white"))
grid.text("\u0394 2099", x=unit(0.68, "npc"), y=unit(0.75, "npc"), gp=gpar(fontsize=25,font=2, col="white"))
grid.text("\u0394 2050", x=unit(0.25, "npc"), y=unit(0.25, "npc"), gp=gpar(fontsize=25,font=2, col="white"))
grid.text("\u0394 2099", x=unit(0.68, "npc"), y=unit(0.25, "npc"), gp=gpar(fontsize=25,font=2, col="white"))
dev.off()
```

```{r climchange, out.width="100%", fig.cap=paste0(sfig_num("climchanges", display="cite"), ": Change in number of hot days and number of days with measureable precipitation from 2010 baseline. Panel (a) presents the 25km x 25km grid cell projections of warming in the United States over this century along the RCP8.5 high emissions scenario, with changes taken from the 2010 baseline. We average downscaled climatic model data from NASA's NEX across the 21 models in the ensemble. The annual number of days above 30$^\\circ{}$C is projected to dramatically increase over the course of this century. Panel (b) also pulls on NASA NEX data to plot the projected change in the number of annual days with measurable precipitation (>1mm) from the 2010 baseline.")}
knitr::include_graphics("./climate_changes.png")
```

```{r}
climchanges <- sfig_num(name = "climchanges", caption = "")
```

<!-- contemporaneous regressions-->

```{r meteo.reg}
main <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r1 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r2 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r3 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r4 <- felm(any_ment ~ cutprcp + cuttmax
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
```

```{r meteo.tab, results='asis'}
invisible(stargazer(main, r1, r2, r3, r4,
          type='latex',
          title = paste0(stab_num("meteo.regs", display="cite"), ':', " Respondent-Day Regressions, Varying Meteorological Controls"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c(gsub('cutprcp', 'PRECIP.DAYS $\\\\in$ ', grep("cutprcp", rownames(main$coef), value = TRUE)), gsub('cuttmax', 'AVG.TMAX $\\\\in$ ', grep("cuttmax", rownames(main$coef), value = TRUE)), "AVG.TRANGE", "AVG.HUMID","AVG.CLOUD","AVG.WIND"),
          add.lines = list(c("City FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes", "Yes", "Yes", "Yes"), 
                           c("State:Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
meteo.regs <- stab_num(name = "meteo.regs", caption = "")
```

```{r par.reg}
brfss[, avg.tmax2 := avg.tmax^2]
brfss[, prcp.days2 := prcp.days^2]

r1 <- felm(any_ment ~ prcp.days + avg.tmax + avg.tmax2 + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
```

```{r par.tab, results='asis'}
invisible(stargazer(r1,
          type='latex',
          title = paste0(stab_num("par.regs", display="cite"), ':', " Respondent-Day Regressions, Parametric Specifications"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c("PRECIP.DAYS","AVG.TMAX","AVG.TMAX²"),
          omit = c("avg.trange","avg.humid","avg.cloud","avg.wind"),
          add.lines = list(c("Meteorological Controls", "Yes", "Yes"),
                           c("City FE", "Yes"),
                           c("Date FE", "Yes"), 
                           c("State:Year FE", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
par.regs <- stab_num(name = "par.regs", caption = "")
```

```{r binday.reg}
binday <- felm(any_ment ~ cutprcp + tmaxneginf + tmax5 + tmax10 + tmax20 + tmax25 + tmax30 + tmax35 + tmaxposinf 
            + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r1 <- felm(any_ment ~ cutprcp + tmaxneginf + tmax5 + tmax10 + tmax20 + tmax25 + tmax30 + tmax35 + tmaxposinf 
            + avg.trange + avg.humid + avg.cloud
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r2 <- felm(any_ment ~ cutprcp + tmaxneginf + tmax5 + tmax10 + tmax20 + tmax25 + tmax30 + tmax35 + tmaxposinf 
            + avg.trange + avg.humid
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r3 <- felm(any_ment ~ cutprcp + tmaxneginf + tmax5 + tmax10 + tmax20 + tmax25 + tmax30 + tmax35 + tmaxposinf 
            + avg.trange
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r4 <- felm(any_ment ~ cutprcp + tmaxneginf + tmax5 + tmax10 + tmax20 + tmax25 + tmax30 + tmax35 + tmaxposinf 
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
```

```{r binday.tab, results='asis'}
invisible(stargazer(binday, r1, r2, r3, r4,
          type='latex',
          title = paste0(stab_num("binday.regs", display="cite"), ':', " Respondent-Day Regressions, Binned Count of Daily Max. Temp."),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c(gsub('cutprcp', 'PRECIP.DAYS $\\\\in$ ', grep("cutprcp", rownames(binday$coef), value = TRUE)), "$\\sum_{r=t-30}^{t} \\mathbf{1}(tmax_{jr} \\in (-\\infty,0])$", "$\\sum_{r=t-30}^{t} \\mathbf{1}(tmax_{jr} \\in (0,5])$", "$\\sum_{r=t-30}^{t} \\mathbf{1}(tmax_{jr} \\in (5,10])$", "$\\sum_{r=t-30}^{t} \\mathbf{1}(tmax_{jr} \\in (15,20])$", "$\\sum_{r=t-30}^{t} \\mathbf{1}(tmax_{jr} \\in (20,25])$", "$\\sum_{r=t-30}^{t} \\mathbf{1}(tmax_{jr} \\in (25,30])$", "$\\sum_{r=t-30}^{t} \\mathbf{1}(tmax_{jr} \\in (30,35])$", "$\\sum_{r=t-30}^{t} \\mathbf{1}(tmax_{jr} \\in (35,\\infty])$", "AVG.TRANGE", "AVG.HUMID","AVG.CLOUD","AVG.WIND"),
          add.lines = list(c("City FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes", "Yes", "Yes", "Yes"), 
                           c("State:Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
binday.regs <- stab_num(name = "binday.regs", caption = "")
```

```{r fx.reg}
brfss[, yrmn := as.factor(paste0(year(time), month(time)))]
brfss[, yrmnwk := as.factor(paste0(year(time), month(time), week(time)))]

r1 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r2 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + state
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r3 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | yrmn + mmsa
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r4 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | yrmnwk + mmsa
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
```

```{r fx.tab, results='asis'}
invisible(stargazer(main, r1, r2, r3, r4,
          type='latex',
          title = paste0(stab_num("fx.regs", display="cite"), ':', " Respondent-Day Regressions, Varying Fixed Effects"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c(gsub('cutprcp', 'PRECIP.DAYS $\\\\in$ ', grep("cutprcp", rownames(main$coef), value = TRUE)), gsub('cuttmax', 'AVG.TMAX $\\\\in$ ', grep("cuttmax", rownames(main$coef), value = TRUE))),
          omit = c("avg.trange","avg.humid","avg.cloud","avg.wind"),
          add.lines = list(c("Meteorological Controls", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("City FE", "Yes", "Yes", "No", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes", "Yes", "No", "No"), 
                           c("State:Year FE", "Yes", "No", "No", "No", "No"),
                           c("State FE", "No", "No", "Yes", "No", "No"),
                           c("Year-Month FE", "No", "No", "No", "Yes", "No"),
                           c("Year-Month-Week FE", "No", "No", "No", "No", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
fx.regs <- stab_num(name = "fx.regs", caption = "")
```

```{r dem.reg}
# demographic fixed effects----
brfss[, cutage := cut(age, breaks=c(-Inf, seq(25, 75, by=10), Inf))] # create factor variable for regresison
brfss[, cutage := relevel(cutage, ref="(45,55]")] # set reference category
brfss[, cutadult := cut(numadult, breaks=c(-Inf, seq(1, 4, by=1), Inf))] # create factor variable
brfss[, cutadult := relevel(cutadult, ref="(1,2]")]
brfss[, cutchild := cut(children, breaks=c(-Inf, seq(1, 5, by=1), Inf))] # create factor variable
brfss[, cutchild := relevel(cutchild, ref="(1,2]")]

# regressions----
r1 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + as.factor(educa) + as.factor(hispanc2) + as.factor(marital) + cutadult + cutchild + as.factor(veteran) + as.factor(bmicat) + time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r2 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + as.factor(educa) + as.factor(hispanc2) + as.factor(marital) + cutadult + cutchild + time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r3 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + as.factor(educa) + as.factor(hispanc2) + time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r4 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r5 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud
          | cutage + as.factor(female) + time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r6 <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + as.factor(educa) + as.factor(hispanc2) + as.factor(marital) + cutadult + cutchild + as.factor(veteran) + as.factor(bmicat) + as.factor(orace2) + time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
```

```{r dem.tab, results='asis'}
invisible(stargazer(r1,r2,r3,r4,r5,r6,
          type='latex',
          title = paste0(stab_num("dem.regs", display="cite"), ':', " Respondent-Day Regressions, Adding Demographic Controls"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c(gsub('cutprcp', 'PRECIP.DAYS $\\\\in$ ', grep("cutprcp", rownames(main$coef), value = TRUE)), gsub('cuttmax', 'AVG.TMAX $\\\\in$ ', grep("cuttmax", rownames(main$coef), value = TRUE))),
          omit=c("avg.trange", "avg.humid", "avg.cloud", "avg.wind"),
          add.lines = list(c("Age FE", "Yes", "Yes","Yes","Yes","Yes","Yes"),
                           c("Gender FE", "Yes", "Yes","Yes","Yes","Yes","Yes"),
                           c("Income FE", "Yes", "Yes","Yes","Yes","No","Yes"),
                           c("Employed FE", "Yes", "Yes","Yes","Yes","No","Yes"),
                           c("Education FE", "Yes", "Yes","Yes","No","No","Yes"),
                           c("Is Hispanic FE", "Yes", "Yes","Yes","No","No","Yes"),
                           c("Marital Status FE", "Yes", "Yes","No","No","No","Yes"),
                           c("Num. Adults in HH FE", "Yes", "Yes","No","No","No","Yes"),
                           c("Num. Children in HH FE", "Yes", "Yes","No","No","No","Yes"),
                           c("Veteran Status FE", "Yes", "No","No","No","No","Yes"),
                           c("BMI Category FE", "Yes", "No","No","No","No","Yes"),
                           c("Race Category FE", "No", "No","No","No","No","Yes"),
                           c("Meteorological Controls", "Yes", "Yes","Yes","Yes","Yes","Yes"),
                           c("City FE", "Yes", "Yes","Yes","Yes","Yes","Yes"),
                           c("Date FE", "Yes", "Yes","Yes","Yes","Yes","Yes"), 
                           c("State:Year FE", "Yes", "Yes","Yes","Yes","Yes","Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
dem.regs <- stab_num(name = "dem.regs", caption = "")
```

```{r city.reg}
city <- brfss[, .(pct_any_ment = mean(any_ment, na.rm=T),
                  avg.trange = unique(avg.trange),
                  avg.cloud = unique(avg.cloud),
                  avg.humid = unique(avg.humid),
                  avg.wind = unique(avg.wind),
                  cutprcp = unique(cutprcp),
                  cuttmax = unique(cuttmax),
                  num = uniqueN(.I)), by=.(mmsa, time, state, stateyr)]

r1 <- felm(pct_any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud
          | time + mmsa + stateyr
          | 0 
          | state
          , data=city,
          weights=city[, num],
          exactDOF=TRUE)
r2 <- felm(pct_any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid
          | time + mmsa + stateyr
          | 0 
          | state
          , data=city,
          weights=city[, num],
          exactDOF=TRUE)
r3 <- felm(pct_any_ment ~ cutprcp + cuttmax + avg.trange
          | time + mmsa + stateyr
          | 0 
          | state
          , data=city,
          weights=city[, num],
          exactDOF=TRUE)
r4 <- felm(pct_any_ment ~ cutprcp + cuttmax
          | time + mmsa + stateyr
          | 0 
          | state
          , data=city,
          weights=city[, num],
          exactDOF=TRUE)
```

```{r city.tab, results='asis'}
invisible(stargazer(r1, r2, r3, r4,
          type='latex',
          title = paste0(stab_num("city.regs", display="cite"), ':', " City-Day Regressions, Varying Meteorological Controls"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = '\\% Reporting Any Mental Health Issue',
          covariate.labels = c(gsub('cutprcp', 'PRECIP.DAYS $\\\\in$ ', grep("cutprcp", rownames(main$coef), value = TRUE)), gsub('cuttmax', 'AVG.TMAX $\\\\in$ ', grep("cuttmax", rownames(main$coef), value = TRUE)), "AVG.TRANGE", "AVG.HUMID","AVG.CLOUD","AVG.WIND"),
          add.lines = list(c("City FE", "Yes", "Yes", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes", "Yes", "Yes"), 
                           c("State:Year FE", "Yes", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Prevalence is reported in percentage points.','Regression is weighted by number of respondents per city-day.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
city.regs <- stab_num(name = "city.regs", caption = "")
```

```{r ext.reg}
# get count of unique respondents per city-day----
ext <- brfss[, .(count = uniqueN(.I),
                 log_count = log(uniqueN(.I)),
                 tmax = unique(tmax),
                 prcp = unique(prcp),
                 humid = unique(humid),
                 cloud = unique(cloud),
                 wind = unique(wind),
                 tmin = unique(tmin),
                 cuttmax = unique(cuttmax),
                 avg.trange = unique(avg.trange),
                 avg.tmax = unique(avg.tmax),
                 cutprcp = unique(cutprcp),
                 prcp.days = unique(prcp.days),
                 avg.humid = unique(avg.humid),
                 avg.cloud = unique(avg.cloud),
                 avg.wind = unique(avg.wind)
                 ), by=.(mmsa, time, stateyr, state, name)]
ext[, avg.tmax2 := avg.tmax^2]

# conduct analysis of 30-day meteorological conditions on extensive margins----
r1 <- felm(log_count ~ prcp.days + avg.tmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=ext,
          exactDOF=TRUE)
r2 <- felm(log_count ~ prcp.days + avg.tmax + avg.tmax2 + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=ext,
          exactDOF=TRUE)

# extensive margin of days a city is in or out of the sample----
setkey(ext, mmsa, time)
setkey(prism, mmsa, time)
ext.full <- merge(prism, subset(ext, select=c('mmsa', 'time', 'count')), all.x=T)
setkey(ext.full, mmsa)
state <- ext[, .(state = unique(state)), by=.(mmsa)]
setkey(state, mmsa)
ext.full <- merge(ext.full, state)
ext.full[, in_sample := !is.na(count)]
ext.full[, stateyr := as.factor(paste0(state, year(time)))]

ext.full[, avg.tmax2 := avg.tmax^2]

## 30-day meteorology
r3 <- felm(in_sample ~ prcp.days + avg.tmax + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=ext.full,
          exactDOF=TRUE)
r4 <- felm(in_sample ~ prcp.days + avg.tmax + avg.tmax2 + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=ext.full,
          exactDOF=TRUE)
```

```{r ext.tab, results='asis'}
invisible(stargazer(r1, r2, r3, r4,
          type='latex',
          title = paste0(stab_num("ext.regs", display="cite"), ':', " Extensive Margin Regressions"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = c('Log \\# Respondents','City in Sample (0/1)'),
          covariate.labels = c("PRECIP.DAYS","AVG.TMAX","AVG.TMAX²"),
          omit = c("avg.trange","avg.humid","avg.cloud","avg.wind"),
          add.lines = list(c("Meteorological Controls", "Yes", "Yes", "Yes", "Yes"),
                           c("City FE", "Yes", "Yes", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes", "Yes", "Yes"), 
                           c("State:Year FE", "Yes", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
ext.regs <- stab_num(name = "ext.regs", caption = "")
```

```{r hetero.reg}
poor <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[income2 <= 22500,],
            exactDOF=TRUE)

rich <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[income2 >= 90000,],
            exactDOF=TRUE)
female <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[female == 1,],
            exactDOF=TRUE)

male <- felm(any_ment ~ cutprcp + cuttmax + avg.trange + avg.humid + avg.cloud + avg.wind
            | time + mmsa + stateyr
            | 0
            | state
            , data=brfss[female == 0,],
            exactDOF=TRUE)
```

```{r hetero.tab, results='asis'}
invisible(stargazer(poor, rich, female, male,
          type='latex',
          title = paste0(stab_num("hetero.regs", display="cite"), ':', " Heterogeneous Effects Regressions"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          column.labels = c('Low Income','High Income','Female','Male'),
          covariate.labels = c(gsub('cutprcp', 'PRECIP.DAYS $\\\\in$ ', grep("cutprcp", rownames(poor$coef), value = TRUE)), gsub('cuttmax', 'AVG.TMAX $\\\\in$ ', grep("cuttmax", rownames(poor$coef), value = TRUE))),
          omit = c("avg.trange","avg.humid","avg.cloud","avg.wind"),
          add.lines = list(c("Meteorological Controls", "Yes", "Yes", "Yes", "Yes"),
                           c("City FE", "Yes", "Yes", "Yes", "Yes"),
                           c("Date FE", "Yes", "Yes", "Yes", "Yes"), 
                           c("State:Year FE", "Yes", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="4pt",
          table.placement = "h"
          ))
```

```{r}
hetero.regs <- stab_num(name = "hetero.regs", caption = "")
```

<!-- long-differences -->

```{r ld.reg}
# long diff----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[year(time) %in% pre.years,
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[year(time) %in% pre.years,
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[year(time) %in% post.years,
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[year(time) %in% post.years,
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN
diff[, tmax.diff2 := tmax.diff^2]

# regressions-----
ld <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample
```

```{r ld.tab, results='asis'}
invisible(stargazer(ld, ld.full, ld.partial,
          type='latex',
          title = paste0(stab_num("ld.regs", display="cite"), ':', " Long Differences Regressions"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Change in \\% Reporting Mental Health Issues',
          column.labels = c('All','No Missing','Missing'),
          notes = c('Standard errors are in parentheses and are clustered on state.','Change is reported in percentage points.'),
          omit = c("prcp.diff","trange.diff","cloud.diff","humid.diff","wind.diff"),
          covariate.labels = c("$\\Delta$ Avg. Max. Temp."),
          add.lines = list(c("Meteorological Controls", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
ld.regs <- stab_num(name = "ld.regs", caption = "")
```

```{r ld.pct.reg}
# long diff----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[year(time) %in% pre.years,
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[year(time) %in% pre.years,
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[year(time) %in% post.years,
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[year(time) %in% post.years,
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN
diff[, tmax.diff2 := tmax.diff^2]
diff[, pct.tmax.diff := 100*((tmax.post - tmax.pre)/tmax.pre)]

ld <- felm(any_ment.diff ~ pct.tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full <- felm(any_ment.diff ~ pct.tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial <- felm(any_ment.diff ~ pct.tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample
```

```{r ld.pct.tab, results='asis'}
invisible(stargazer(ld, ld.full, ld.partial,
          type='latex',
          title = paste0(stab_num("ld.pct.regs", display="cite"), ':', " Long Differences Regressions, \\% Change in Max. Temp."),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Change in \\% Reporting Mental Health Issues',
          column.labels = c('All','No Missing','Missing'),
          notes = c('Standard errors are in parentheses and are clustered on state.','Change is reported in percentage points.'),
          omit = c("prcp.diff","trange.diff","cloud.diff","humid.diff","wind.diff"),
          covariate.labels = c("$\\% \\Delta$ Avg. Max. Temp."),
          add.lines = list(c("Meteorological Controls", "Yes", "Yes", "Yes")),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
ld.pct.regs <- stab_num(name = "ld.pct.regs", caption = "")
```

```{r ld.dem.reg}
# long diff----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[year(time) %in% pre.years,
               .(any_ment = mean(any_ment, na.rm=T),
                 age = mean(age, na.rm=T),
                 hispanc2 = mean(hispanc2, na.rm=T),
                 highschool = mean(educa==3|educa==4, na.rm=T),
                 college = mean(educa==5|educa==6, na.rm=T),
                 income2 = mean(income2, na.rm=T),
                 employ = mean(employ, na.rm=T),
                 female = mean(female, na.rm=T),
                 bmi = mean(bmi, na.rm=T),
                 married = mean(marital==1, na.rm=T),
                 divorced = mean(marital==2, na.rm=T),
                 white = mean(orace2==1, na.rm=T),
                 black = mean(orace2==2, na.rm=T),
                 numadult = mean(numadult, na.rm=T),
                 children = mean(children, na.rm=T),
                 veteran = mean(veteran, na.rm=T),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon)]
setnames(r.pre, c("any_ment","age","hispanc2","highschool","college","income2","employ","female","bmi","married","divorced","white","black","numadult","children","veteran","years", "subs"), paste0(c("any_ment","age","hispanc2","highschool","college","income2","employ","female","bmi","married","divorced","white","black","numadult","children","veteran","years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[year(time) %in% pre.years,
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid", "wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid", "wind"), paste0(c("tmax","trange","prcp","cloud","humid", "wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[year(time) %in% post.years,
               .(any_ment = mean(any_ment, na.rm=T),
                 age = mean(age, na.rm=T),
                 hispanc2 = mean(hispanc2, na.rm=T),
                 highschool = mean(educa==3|educa==4, na.rm=T),
                 college = mean(educa==5|educa==6, na.rm=T),
                 income2 = mean(income2, na.rm=T),
                 employ = mean(employ, na.rm=T),
                 female = mean(female, na.rm=T),
                 bmi = mean(bmi, na.rm=T),
                 married = mean(marital=='married', na.rm=T),
                 divorced = mean(marital=='divorced', na.rm=T),
                 white = mean(orace2==1, na.rm=T),
                 black = mean(orace2==2, na.rm=T),
                 numadult = mean(numadult, na.rm=T),
                 children = mean(children, na.rm=T),
                 veteran = mean(veteran, na.rm=T),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon)]
setnames(r.post, c("any_ment","age","hispanc2","highschool","college","income2","employ","female","bmi","married","divorced","white","black","numadult","children","veteran","years", "subs"), paste0(c("any_ment","age","hispanc2","highschool","college","income2","employ","female","bmi","married","divorced","white","black","numadult","children","veteran","years", "subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[year(time) %in% post.years,
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid", "wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid", "wind"), paste0(c("tmax","trange","prcp","cloud","humid", "wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

# regressions----
ld1 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
           + age.diff + female.diff + income2.diff + employ.diff + highschool.diff + college.diff + hispanc2.diff + married.diff + divorced.diff + numadult.diff + children.diff + veteran.diff + bmi.diff
           | 0 
           | 0 
           | state, 
           data=diff, 
           weights=diff[, subs])
ld2 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
           + age.diff + female.diff + income2.diff + employ.diff + highschool.diff + college.diff + hispanc2.diff + married.diff + divorced.diff + numadult.diff + children.diff
           | 0 
           | 0 
           | state, 
           data=diff, 
           weights=diff[, subs])
ld3 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
           + age.diff + female.diff + income2.diff + employ.diff + highschool.diff + college.diff + hispanc2.diff
           | 0 
           | 0 
           | state, 
           data=diff, 
           weights=diff[, subs])
ld4 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
           + age.diff + female.diff + income2.diff + employ.diff
           | 0 
           | 0 
           | state, 
           data=diff, 
           weights=diff[, subs])
ld5 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
           + age.diff + female.diff
           | 0 
           | 0 
           | state, 
           data=diff, 
           weights=diff[, subs])
ld6 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
           + age.diff + female.diff + income2.diff + employ.diff + highschool.diff + college.diff + hispanc2.diff + married.diff + divorced.diff + numadult.diff + children.diff + veteran.diff + bmi.diff + white.diff + black.diff
           | 0 
           | 0 
           | state, 
           data=diff, 
           weights=diff[, subs])
```

```{r ld.dem.tab, results='asis'}
invisible(stargazer(ld1,ld2,ld3,ld4,ld5,ld6,
          type='latex',
          title = paste0(stab_num("ld.dem.regs", display="cite"), ':', " All Cities Long Differences Regressions, Demographic Controls"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Change in \\% Reporting Mental Health Issues',
          notes = c('Standard errors are in parentheses and are clustered on state.','Change is reported in percentage points.'),
          omit = c("prcp.diff","trange.diff","cloud.diff","humid.diff","wind.diff", grep('age', rownames(ld6$coef), value=TRUE), grep('female', rownames(ld6$coef), value=TRUE), grep('income', rownames(ld6$coef), value=TRUE), grep('employ', rownames(ld6$coef), value=TRUE), grep('highschool', rownames(ld6$coef), value=TRUE), grep('college', rownames(ld6$coef), value=TRUE), grep('hispanc', rownames(ld6$coef), value=TRUE), grep('married', rownames(ld6$coef), value=TRUE), grep('divorced', rownames(ld6$coef), value=TRUE), grep('numadult', rownames(ld6$coef), value=TRUE), grep('children', rownames(ld6$coef), value=TRUE), grep('veteran', rownames(ld6$coef), value=TRUE), grep('bmi', rownames(ld6$coef), value=TRUE), grep('white', rownames(ld6$coef), value=TRUE), grep('black', rownames(ld6$coef), value=TRUE)),
          covariate.labels = c("$\\Delta$ Avg. Max. Temp."),
          add.lines = list(c("$\\Delta$ Mean Age", "Yes", "Yes","Yes","Yes","Yes","Yes"),
                           c("$\\Delta$ \\% Female", "Yes", "Yes","Yes","Yes","Yes","Yes"),
                           c("$\\Delta$ Mean Income", "Yes", "Yes","Yes","Yes","No","Yes"),
                           c("$\\Delta$ \\% Employed", "Yes", "Yes","Yes","Yes","No","Yes"),
                           c("$\\Delta$ \\% High School Education", "Yes", "Yes","Yes","No","No","Yes"),
                           c("$\\Delta$ \\% College Education", "Yes", "Yes","Yes","No","No","Yes"),
                           c("$\\Delta$ \\% Hispanic", "Yes", "Yes","Yes","No","No","Yes"),
                           c("$\\Delta$ \\% Married", "Yes", "Yes","No","No","No","Yes"),
                           c("$\\Delta$ \\% Divorced", "Yes", "Yes","No","No","No","Yes"),
                           c("$\\Delta$ Mean Num. Adults in HH", "Yes", "Yes","No","No","No","Yes"),
                           c("$\\Delta$ Mean Num. Childen in HH", "Yes", "Yes","No","No","No","Yes"),
                           c("$\\Delta$ \\% Veteran", "Yes", "No","No","No","No","Yes"),
                           c("$\\Delta$ Mean BMI", "Yes", "No","No","No","No","Yes"),
                           c("$\\Delta$ \\% White", "No", "No","No","No","No","Yes"),
                           c("$\\Delta$ \\% Black", "No", "No","No","No","No","Yes"),
                           c("Meteorological Controls", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
ld.dem.regs <- stab_num(name = "ld.dem.regs", caption = "")
```

```{r ld.census.dem.reg}
# long diff----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[year(time) %in% pre.years,
               .(any_ment = mean(any_ment, na.rm=T),
                 age = mean(age, na.rm=T),
                 hispanc2 = mean(hispanc2, na.rm=T),
                 highschool = mean(educa==3|educa==4, na.rm=T),
                 college = mean(educa==5|educa==6, na.rm=T),
                 income2 = mean(income2, na.rm=T),
                 employ = mean(employ, na.rm=T),
                 female = mean(female, na.rm=T),
                 bmi = mean(bmi, na.rm=T),
                 married = mean(marital==1, na.rm=T),
                 divorced = mean(marital==2, na.rm=T),
                 white = mean(orace2==1, na.rm=T),
                 black = mean(orace2==2, na.rm=T),
                 numadult = mean(numadult, na.rm=T),
                 children = mean(children, na.rm=T),
                 veteran = mean(veteran, na.rm=T),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon)]
setnames(r.pre, c("any_ment","age","hispanc2","highschool","college","income2","employ","female","bmi","married","divorced","white","black","numadult","children","veteran","years", "subs"), paste0(c("any_ment","age","hispanc2","highschool","college","income2","employ","female","bmi","married","divorced","white","black","numadult","children","veteran","years", "subs"), '.pre'))
setkey(r.pre, mmsa, state)

## census pre-period averages
s.pre <- seer[year %in% pre.years,
               .(pct.female = mean(pct.female, na.rm=T),
                 pct.white = mean(pct.white, na.rm=T),
                 pct.black = mean(pct.black, na.rm=T),
                 mean.pop = mean(totpop, na.rm=T),
                 median.age = mean(median.age, na.rm=T)),
               by=.(mmsa, state)]
setnames(s.pre, c("pct.female","pct.white","pct.black","mean.pop","median.age"), paste0(c("pct.female","pct.white","pct.black","mean.pop","median.age"), '.pre'))
setkey(s.pre, mmsa, state)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[year(time) %in% pre.years,
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid", "wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid", "wind"), paste0(c("tmax","trange","prcp","cloud","humid", "wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, s.pre)
setkey(pre, mmsa)
pre <- merge(pre, c.pre)
rm(r.pre, c.pre, s.pre)

## respondent post-period averages
r.post <- brfss[year(time) %in% post.years,
               .(any_ment = mean(any_ment, na.rm=T),
                 age = mean(age, na.rm=T),
                 hispanc2 = mean(hispanc2, na.rm=T),
                 highschool = mean(educa==3|educa==4, na.rm=T),
                 college = mean(educa==5|educa==6, na.rm=T),
                 income2 = mean(income2, na.rm=T),
                 employ = mean(employ, na.rm=T),
                 female = mean(female, na.rm=T),
                 bmi = mean(bmi, na.rm=T),
                 married = mean(marital=='married', na.rm=T),
                 divorced = mean(marital=='divorced', na.rm=T),
                 white = mean(orace2==1, na.rm=T),
                 black = mean(orace2==2, na.rm=T),
                 numadult = mean(numadult, na.rm=T),
                 children = mean(children, na.rm=T),
                 veteran = mean(veteran, na.rm=T),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon)]
setnames(r.post, c("any_ment","age","hispanc2","highschool","college","income2","employ","female","bmi","married","divorced","white","black","numadult","children","veteran","years", "subs"), paste0(c("any_ment","age","hispanc2","highschool","college","income2","employ","female","bmi","married","divorced","white","black","numadult","children","veteran","years", "subs"), '.post'))
setkey(r.post, mmsa, state)

## census post-period averages
s.post <- seer[year %in% post.years,
               .(pct.female = mean(pct.female, na.rm=T),
                 pct.white = mean(pct.white, na.rm=T),
                 pct.black = mean(pct.black, na.rm=T),
                 mean.pop = mean(totpop, na.rm=T),
                 median.age = mean(median.age, na.rm=T)),
               by=.(mmsa, state)]
setnames(s.post, c("pct.female","pct.white","pct.black","mean.pop","median.age"), paste0(c("pct.female","pct.white","pct.black","mean.pop","median.age"), '.post'))
setkey(s.post, mmsa, state)

## city post-period averages
c.post <- prism[year(time) %in% post.years,
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid", "wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid", "wind"), paste0(c("tmax","trange","prcp","cloud","humid", "wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, s.post)
setkey(post, mmsa)
post <- merge(post, c.post)
rm(r.post, s.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

# regressions----
ld1 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
            + median.age.diff + pct.female.diff + mean.pop.diff + pct.white.diff + pct.black.diff
            | 0 
            | 0 
            | state, 
            data=diff, 
            weights=diff[, subs])
ld2 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
            + median.age.diff + pct.female.diff + mean.pop.diff
            | 0 
            | 0 
            | state, 
            data=diff, 
            weights=diff[, subs])
ld3 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
            + median.age.diff + pct.female.diff
            | 0 
            | 0 
            | state, 
            data=diff, 
            weights=diff[, subs])
ld4 <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff
            + median.age.diff
            | 0 
            | 0 
            | state, 
            data=diff, 
            weights=diff[, subs])
```

```{r ld.census.dem.tab, results='asis'}
invisible(stargazer(ld1,ld2,ld3,ld4,
          type='latex',
          title = paste0(stab_num("ld.census.dem.regs", display="cite"), ':', " All Cities Long Differences, Census Demographic Controls"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Change in \\% Reporting Mental Health Issues',
          notes = c('Standard errors are in parentheses and are clustered on state.','Change is reported in percentage points.'),
          omit = c("prcp.diff","trange.diff","cloud.diff","humid.diff","wind.diff", grep('age', rownames(ld1$coef), value=TRUE), grep('female', rownames(ld1$coef), value=TRUE), grep('pop', rownames(ld1$coef), value=TRUE), grep('white', rownames(ld1$coef), value=TRUE), grep('black', rownames(ld1$coef), value=TRUE)),
          covariate.labels = c("$\\Delta$ Avg. Max. Temp."),
          add.lines = list(c("$\\Delta$ Median Age", "Yes", "Yes","Yes","Yes"),
                           c("$\\Delta$ \\% Female", "Yes", "Yes","Yes","No"),
                           c("$\\Delta$ Population", "Yes", "Yes","No","No"),
                           c("$\\Delta$ \\% White", "Yes", "No","No","No"),
                           c("$\\Delta$ \\% Black", "Yes", "No","No","No"),
                           c("Meteorological Controls", "Yes", "Yes", "Yes", "Yes")),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
ld.census.dem.regs <- stab_num(name = "ld.census.dem.regs", caption = "")
```

```{r ld.seas.reg}
# long-diff spring----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[(month(time) %in% c(3,4,5)) & (year(time) %in% pre.years),
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[(month(time) %in% c(3,4,5)) & (year(time) %in% pre.years),
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[(month(time) %in% c(3,4,5)) & (year(time) %in% post.years),
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[(month(time) %in% c(3,4,5)) & (year(time) %in% post.years),
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld.spr <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full.spr <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial.spr <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample

# long-diff summer----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[(month(time) %in% 6:8) & (year(time) %in% pre.years),
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[(month(time) %in% 6:8) & (year(time) %in% pre.years),
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[(month(time) %in% 6:8) & (year(time) %in% post.years),
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[(month(time) %in% 6:8) & (year(time) %in% post.years),
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld.sum <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full.sum <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial.sum <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample

# long-diff fall----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[(month(time) %in% c(9,10,11)) & (year(time) %in% pre.years),
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[(month(time) %in% c(9,10,11)) & (year(time) %in% pre.years),
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[(month(time) %in% c(9,10,11)) & (year(time) %in% post.years),
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[(month(time) %in% c(9,10,11)) & (year(time) %in% post.years),
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld.fall <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full.fall <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial.fall <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample

# long-diff winter----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[(month(time) %in% c(12,1,2)) & (year(time) %in% pre.years),
               .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("any_ment")]
setnames(r.pre, c("any_ment", "years", "subs"), paste0(c("any_ment", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[(month(time) %in% c(12,1,2)) & (year(time) %in% pre.years),
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[(month(time) %in% c(12,1,2)) & (year(time) %in% post.years),
                .(any_ment = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("any_ment")]
setnames(r.post, c("any_ment","years","subs"), paste0(c("any_ment","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[(month(time) %in% c(12,1,2)) & (year(time) %in% post.years),
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN

ld.win <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])
ld.full.win <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years==10,], weights=diff[years==10, subs]) # cities with all 10 years in sample
ld.partial.win <- felm(any_ment.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff[years<10,], weights=diff[years<10, subs]) # cities will less than 10 years in sample
```

```{r ld.seas.tab, results='asis'}
invisible(stargazer(ld.spr, ld.sum, ld.fall, ld.win,
          type='latex',
          title = paste0(stab_num("ld.seas.regs", display="cite"), ':', " Long Differences Regressions by Seasons"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Change in \\% Reporting Mental Health Issues',
          column.labels = c('Spring','Summer','Fall', 'Winter'),
          notes = c('Standard errors are in parentheses and are clustered on state.','Change is reported in percentage points.'),
          omit = c("prcp.diff","trange.diff","cloud.diff","humid.diff","wind.diff"),
          covariate.labels = c("$\\Delta$ Avg. Max. Temp."),
          add.lines = list(c("Meteorological Controls", "Yes", "Yes", "Yes", "Yes")),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
ld.seas.regs <- stab_num(name = "ld.seas.regs", caption = "")
```

<!-- dnd -->

```{r dnd.reg}
# diff-n-diff----

## set up treatment variables
brfss[, treat := (state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")] 

## set pre-Katrina as pre and post-Katrina as post
brfss[, post := time > "2005-08-24"]

## interaction
brfss[, postXtreat := post*treat]

## create exposure variable
brfss[, katrina.exposure := ((state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")) & time > "2005-08-24"]

main <- felm(any_ment ~ post + treat + postXtreat
             | 0 
             | 0 
             | state, 
             data=brfss, 
             exactDOF=TRUE)
r1 <- felm(any_ment ~ katrina.exposure 
           | mmsa + time 
           | 0 
           | state, 
           data=brfss, 
           exactDOF=TRUE)
```

```{r dnd.tab, results='asis'}
invisible(stargazer(main, r1,
          type='latex',
          title = paste0(stab_num("dnd.regs", display="cite"), ':', " Katrina Difference-in-Differences Regressions"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c("Post Period", "Katrina Treatment", "Post X Katrina", "Katrina Exposure"),
          add.lines = list(c("City FE", "No", "Yes"),
                           c("Date FE", "No", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
dnd.regs <- stab_num(name = "dnd.regs", caption = "")
```

```{r dnd.dem.reg}
# diff-n-diff, dem controls----

## set up treatment variables
brfss[, treat := (state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")] 

## set pre-Katrina as pre and post-Katrina as post
brfss[, post := time > "2005-08-24"]

## interaction
brfss[, postXtreat := post*treat]

## create exposure variable
brfss[, katrina.exposure := ((state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")) & time > "2005-08-24"]

## demographic fixed effects
brfss[, cutage := cut(age, breaks=c(-Inf, seq(25, 75, by=10), Inf))] # create factor variable for regresison
brfss[, cutage := relevel(cutage, ref="(45,55]")] # set reference category
brfss[, cutadult := cut(numadult, breaks=c(-Inf, seq(1, 4, by=1), Inf))] # create factor variable
brfss[, cutadult := relevel(cutadult, ref="(1,2]")]
brfss[, cutchild := cut(children, breaks=c(-Inf, seq(1, 5, by=1), Inf))] # create factor variable
brfss[, cutchild := relevel(cutchild, ref="(1,2]")]

# regressions----
r1 <- felm(any_ment ~ post + treat + postXtreat
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + as.factor(educa) + as.factor(hispanc2) + as.factor(marital) + cutadult + cutchild + as.factor(veteran) + as.factor(bmicat)
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r2 <- felm(any_ment ~ post + treat + postXtreat
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + as.factor(educa) + as.factor(hispanc2) + as.factor(marital) + cutadult + cutchild
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r3 <- felm(any_ment ~ post + treat + postXtreat
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + as.factor(educa) + as.factor(hispanc2)
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r4 <- felm(any_ment ~ post + treat + postXtreat
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ)
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r5 <- felm(any_ment ~ post + treat + postXtreat
          | cutage + as.factor(female)
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
r6 <- felm(any_ment ~ post + treat + postXtreat
          | cutage + as.factor(female) + as.factor(income2) + as.factor(employ) + as.factor(educa) + as.factor(hispanc2) + as.factor(marital) + cutadult + cutchild + as.factor(veteran) + as.factor(bmicat) + as.factor(orace2)
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)
```

```{r dnd.dem.tab, results='asis'}
invisible(stargazer(r1,r2,r3,r4,r5,r6,
          type='latex',
          title = paste0(stab_num("dnd.dem.regs", display="cite"), ':', " Katrina Difference-in-Differences Regressions, Demographic Controls"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c("Post Period", "Katrina Treatment", "Post X Katrina"),
          add.lines = list(c("Age FE", "Yes", "Yes","Yes","Yes","Yes","Yes"),
                           c("Gender FE", "Yes", "Yes","Yes","Yes","Yes","Yes"),
                           c("Income FE", "Yes", "Yes","Yes","Yes","No","Yes"),
                           c("Employed FE", "Yes", "Yes","Yes","Yes","No","Yes"),
                           c("Education FE", "Yes", "Yes","Yes","No","No","Yes"),
                           c("Is Hispanic FE", "Yes", "Yes","Yes","No","No","Yes"),
                           c("Marital Status FE", "Yes", "Yes","No","No","No","Yes"),
                           c("Num. Adults in HH FE", "Yes", "Yes","No","No","No","Yes"),
                           c("Num. Children in HH FE", "Yes", "Yes","No","No","No","Yes"),
                           c("Veteran Status FE", "Yes", "No","No","No","No","Yes"),
                           c("BMI Category FE", "Yes", "No","No","No","No","Yes"),
                           c("Race Category FE", "No", "No","No","No","No","Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
dnd.dem.regs <- stab_num(name = "dnd.dem.regs", caption = "")
```

```{r dnd.city.reg}
city <- brfss[, .(pct_any_ment = mean(any_ment, na.rm=T),
                  avg.trange = unique(avg.trange),
                  avg.cloud = unique(avg.cloud),
                  avg.humid = unique(avg.humid),
                  avg.wind = unique(avg.wind),
                  cutprcp = unique(cutprcp),
                  cuttmax = unique(cuttmax),
                  num = uniqueN(.I)), by=.(mmsa, name, time, state, stateyr)]

## set up treatment variables
city[, treat := (state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")] 

## set pre-Katrina as pre and post-Katrina as post
city[, post := time > "2005-08-24"]

## interaction
city[, postXtreat := post*treat]

## create exposure variable
city[, katrina.exposure := ((state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")) & time > "2005-08-24"]

r1 <- felm(pct_any_ment ~ post + treat + postXtreat
             | 0 
             | 0 
             | state, 
             data=city,
             weights=city[, num],
             exactDOF=TRUE)
r2 <- felm(pct_any_ment ~ katrina.exposure 
           | mmsa + time 
           | 0 
           | state, 
           data=city,
           weights=city[, num],
           exactDOF=TRUE)
```

```{r dnd.city.tab, results='asis'}
invisible(stargazer(r1, r2,
          type='latex',
          title = paste0(stab_num("dnd.city.regs", display="cite"), ':', " City-Day Difference-in-Differences Regressions"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = '\\% Reporting Any Mental Health Issue',
          covariate.labels = c("Post Period", "Katrina Treatment", "Post X Katrina", "Katrina Exposure"),
          add.lines = list(c("City FE", "No", "Yes"),
                           c("Date FE", "No", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Prevalence is reported in percentage points.','Regression is weighted by number of respondents per city-day.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
dnd.city.regs <- stab_num(name = "dnd.city.regs", caption = "")
```

```{r dnd.alt.reg}
# diff-n-diff using only neighboring cities (without disaster declarations) as controls----

## set up treatment variables
brfss[, treat := (state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")] 

## set pre-Katrina as pre and post-Katrina as post
brfss[, post := time > "2005-08-24"]

## interaction
brfss[, postXtreat := post*treat]

## only look at neighboring states as controls (geographic proximity)
brfss[state %like% "FL|AL|AR|TX|GA|TN" & treat!=TRUE , treat := FALSE]
brfss[!(state %like% "FL|AL|AR|TX|LA|MS|GA|TN"), treat := NA]

## create exposure variable only for affected and neighboring states
brfss[, katrina.exposure := ((state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")) & time > "2005-08-24"]
brfss[is.na(treat), katrina.exposure := NA]

neighbor1 <- felm(any_ment ~ post + treat + postXtreat
             | 0 
             | 0 
             | state, 
             data=brfss, 
             exactDOF=TRUE)
neighbor2 <- felm(any_ment ~ katrina.exposure 
           | mmsa + time 
           | 0 
           | state, 
           data=brfss, 
           exactDOF=TRUE)
```

```{r dnd.alt.tab, results='asis'}
invisible(stargazer(neighbor1, neighbor2,
          type='latex',
          title = paste0(stab_num("dnd.alt.regs", display="cite"), ':', " Katrina Difference-in-Differences Regressions, Neighbors Only"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c("Post Period", "Katrina Treatment", "Post X Katrina", "Katrina Exposure"),
          add.lines = list(c("City FE", "No", "Yes"),
                           c("Date FE", "No", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
dnd.alt.regs <- stab_num(name = "dnd.alt.regs", caption = "")
```

```{r dnd.bal.reg}
# examine diff-n-diff using only the cities that appear in both pre and post periods----

## set up treatment variables
brfss[, treat := (state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")]

## set pre-Katrina as pre and post-Katrina as post
brfss[, post := time > "2005-08-24"]

## interaction
brfss[, postXtreat := post*treat]

## keep only balanced panel of cities in treat/control
pre.cities <- brfss[post == FALSE, unique(name)]
post.cities <- brfss[post == TRUE, unique(name)]
brfss[!(name %in% pre.cities & name %in% post.cities), treat := NA]

## create exposure variable, limit to only balanced panel of cities
brfss[name %in% pre.cities & name %in% post.cities, katrina.exposure := ((state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")) & time > "2005-08-24"]

balanced1 <- felm(any_ment ~ post + treat + postXtreat
                  | 0
                  | 0
                  | state,
                  data=brfss,
                  exactDOF=TRUE)
balanced2 <- felm(any_ment ~ katrina.exposure
                  | mmsa + time
                  | 0
                  | state,
                  data=brfss,
                  exactDOF=TRUE)
```

```{r dnd.bal.tab, results='asis'}
invisible(stargazer(balanced1, balanced2,
          type='latex',
          title = paste0(stab_num("dnd.bal.regs", display="cite"), ':', " Katrina Difference-in-Differences Regressions, Balanced Cities"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = 'Any Mental Health Issue',
          covariate.labels = c("Post Period", "Katrina Treatment", "Post X Katrina", "Katrina Exposure"),
          add.lines = list(c("City FE", "No", "Yes"),
                           c("Date FE", "No", "Yes")),
          notes = c('Standard errors are in parentheses and are clustered on state.','Probability is reported in percentage points.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="20pt",
          table.placement = "h"
          ))
```

```{r}
dnd.bal.regs <- stab_num(name = "dnd.bal.regs", caption = "")
```

<!-- negative binomial -->

```{r rep.reg}
# contemporaneous ols----
brfss[, avg.tmax2 := avg.tmax^2]
rep.quad <- felm(menthlth ~ prcp.days + avg.tmax + avg.tmax2 + avg.trange + avg.humid + avg.cloud + avg.wind
          | time + mmsa + stateyr
          | 0 
          | state
          , data=brfss,
          exactDOF=TRUE)

# contemporaneous negative binomial----
library(mfx)
# # because of size of matrix, can't run full ols model with negative binomial
# # use a other non-parametric controls
# brfss[, avg.tmax2 := avg.tmax^2]
# brfss[, yrmon := as.factor(as.yearmon(time))]
# brfss[, dow := as.factor(weekdays(time))]
# 
# # glm.nb can't do easy clusters, use mfx, but mfx doesn't work with stargazer.
# 
# negbin2.cl <- negbinirr(formula = menthlth ~
#                         prcp.days + avg.tmax + avg.tmax2
#                         + avg.trange + avg.humid + avg.cloud + avg.wind
#                         + dow + yrmon + mmsa,
#                         data=brfss, clustervar1="state")
# negbin2 <- glm.nb(formula = menthlth ~
#                         prcp.days + avg.tmax + avg.tmax2
#                         + avg.trange + avg.humid + avg.cloud + avg.wind
#                         + dow + yrmon + mmsa,
#                         data=brfss)
# 
# rm(list=setdiff(ls(), c('negbin2.cl','negbin2')))
# save.image("./negative_binomial.Rdata", compress=TRUE)

load("./negative_binomial.Rdata")

rep.nb <- negbin2
nbcl.se <- negbin2.cl$irr[, 2]
nbcl.z <- negbin2.cl$irr[, 3]
nbcl.p <- negbin2.cl$irr[, 4]

rm(negbin2.cl, negbin2)
invisible(gc())

# long diff----
pre.years <- 2002:2006
post.years <- 2007:2011

## respondent pre-period averages
r.pre <- brfss[year(time) %in% pre.years,
               .(menthlth = unlist(lapply(.SD, mean, na.rm=T)),
                 years = uniqueN(year(time)),
                 subs = uniqueN(.I)),
               by=.(mmsa, state, lat, lon),
               .SDcols = c("menthlth")]
setnames(r.pre, c("menthlth", "years", "subs"), paste0(c("menthlth", "years", "subs"), '.pre'))
setkey(r.pre, mmsa)

## city pre-period averages
prism[, trange := tmax - tmin]
c.pre <- prism[year(time) %in% pre.years,
               lapply(.SD, mean, na.rm=T),
               by=.(mmsa),
               .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.pre, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.pre'))
setkey(c.pre, mmsa)

## merge pre
pre <- merge(r.pre, c.pre)
rm(r.pre, c.pre)

## respondent post-period averages
r.post <- brfss[year(time) %in% post.years,
                .(menthlth = unlist(lapply(.SD, mean, na.rm=T)),
                  years = uniqueN(year(time)),
                  subs = uniqueN(.I)),
                by=.(mmsa, state, lat, lon),
                .SDcols = c("menthlth")]
setnames(r.post, c("menthlth","years","subs"), paste0(c("menthlth","years","subs"), '.post'))
setkey(r.post, mmsa)

## city post-period averages
c.post <- prism[year(time) %in% post.years,
                lapply(.SD, mean, na.rm=T),
                by=.(mmsa),
                .SDcols = c("tmax","trange","prcp","cloud","humid","wind")]
setnames(c.post, c("tmax","trange","prcp","cloud","humid","wind"), paste0(c("tmax","trange","prcp","cloud","humid","wind"), '.post'))
setkey(c.post, mmsa)

## merge post
post <- merge(r.post, c.post)
rm(r.post, c.post)

## merge
setkey(pre, mmsa, state, lat, lon)
setkey(post, mmsa, state, lat, lon)
diff <- merge(pre, post) # merge together

## calculate total number of years that enter sample
diff[, years := years.pre + years.post]
diff[, subs := subs.pre + subs.post]

## list of variables to difference
v <- gsub('.pre', '', grep('.pre', names(pre), value=TRUE))
invisible(
    foreach::foreach (l=v) %do% {
        diff <- diff[, paste0(l, '.diff') := get(paste0(l, '.post')) - get(paste0(l, '.pre'))]
    } # generate all differences variables
)

diff[, c('years.diff','subs.diff') := NULL]
diff <- diff[!is.nan(tmax.diff),] # remove obs with NaN
diff[, tmax.diff2 := tmax.diff^2]

rep.ld <- felm(menthlth.diff ~ tmax.diff + prcp.diff + trange.diff + cloud.diff + humid.diff + wind.diff | 0 | 0 | state, data=diff, weights=diff[, subs])

# katrina diff-n-diff----

## set up treatment variables
brfss[, treat := (state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")] 

## set pre-Katrina as pre and post-Katrina as post
brfss[, post := time > "2005-08-24"]

## interaction
brfss[, postXtreat := post*treat]

## create exposure variable
brfss[, katrina.exposure := ((state %like% "LA|MS") | (name %like% "Birmingham AL|Mobile AL|Tuscaloosa AL|Miami FL|Key West FL|Naples FL|Panama City FL|Pensacola FL")) & time > "2005-08-24"]

rep.diffndiff <- felm(menthlth ~ post + treat + postXtreat
                         | 0 
                         | 0 
                         | state, 
                         data=brfss, 
                         exactDOF=TRUE)
```

```{r rep.tab, results='asis'}
invisible(stargazer(rep.quad, rep.nb, rep.ld, rep.diffndiff,
          se = list(NULL, nbcl.se, NULL, NULL),
          p = list(NULL, nbcl.p, NULL, NULL),
          t = list(NULL, nbcl.z, NULL, NULL),
          type='latex',
          title = paste0(stab_num("rep.regs", display="cite"), ':', " Regressions with Count of Poor Mental Health Days"),
          model.numbers = TRUE,
          dep.var.labels.include = TRUE,
          dep.var.labels = c('Mental Health Days','$\\Delta$ Mental Health Days','Mental Health Days'),
          covariate.labels = c("PRECIP.DAYS","AVG.TMAX","AVG.TMAX²","$\\Delta$ Avg. Max. Temp.","Post Period", "Katrina Treatment", "Post X Katrina"),
          omit = c("avg.trange","avg.humid","avg.cloud","avg.wind","prcp.diff","trange.diff","humid.diff","cloud.diff","wind.diff", glob2rx("mmsa*"),glob2rx("dow*"),glob2rx("yrmon*")),
          add.lines = list(c("Meteorological Controls", "Yes", "Yes", "Yes", "No"),
                           c("City FE", "Yes", "Yes", "No", "No"),
                           c("Date FE", "Yes", "No", "No", "No"), 
                           c("State:Year FE", "Yes", "No", "No", "No"),
                           c("Calendar Month FE", "No", "Yes", "No", "No"),
                           c("Day of Week FE", "No", "Yes", "No", "No")),
          notes = c('Standard errors are in parentheses and are clustered on state.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size="footnotesize",
          column.sep.width="0pt",
          table.placement = "h"
          ))
```

```{r}
rep.regs <- stab_num(name = "rep.regs", caption = "")
```
